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Abstract 

We examine the dynamics of the translation stage of cellular protein production, in which ribo- 
somes move uni-directionally along mRNA strands building an amino acid chain as they go. We 
describe the system using a timed event graph - a class of Petri net useful for studying discrete 
events which take a finite time. We use max-plus algebra to describe a deterministic version of the 
model, calculating the protein production rate and density of ribosomes on the mRNA. We find exact 
agreement between these analytical results and numerical simulations of the deterministic case. 

1 Introduction 

Messenger RNA (mRNA) translation is one of the steps in protein production in cells pQ. mRNAs are 
single strands of nucleotides which are transcribed from the DNA. The sequence of nucleotides, grouped 
in triplets called codons, holds the code for a specific chain of amino acids that makes up a protein. 
Translation is performed by molecular machines called ribosomes, which scan along the mRNA adding 
amino acids to a growing chain which will become the protein. The rates at which different proteins are 
produced are crucial in determining how a cell grows and functions. 

Various statistical models have been used to describe and understand the translation process. In this 
paper we propose a new model which can be analysed rather completely using algebra on the max-plus 
semi-ring. The most convenient way to present the model is to write it as a deterministic Petri net [2]. 
This can then be analysed completely as a linear dynamical system when written in terms of the max-plus 
algebra [3 . The model gives a more realistic representation of the movement of ribosomes than that of 
many previous studies, and explicitly considers the exact genetic coding sequence as well as the fact that 
different codons are translated at different rates. We perform numerical simulations of the Petri net, and 
find the results to exactly match the behaviour predicted by the max-plus algebra. 

Although the general sequence of events in translation is understood, there are many questions which 
remain unanswered. Translation begins when a ribosome binds to the end of an mRNA strand - a 
process known as initiation. In this paper we focus on the second stage of translation, which is known 
as elongation. Here the ribosome moves uni-directionally along the mRNA, pausing at each codon to 
recruit an amino acid, which is then added to the growing chain. Amino acids are transported within the 
cytoplasm by carrier molecules called transfer RNAs (tRNAs). Different species of tRNA carry different 
amino acids, and different codons correspond to different tRNAs. The different tRNAs appear in different 
concentrations, so the time which a ribosome waits for the required tRNA to arrive differs from codon to 
codon [H [5] . Since it is often the case that several ribosomes are bound to the mRNA at the same time, 
the distribution of waiting times can result in traffic jams if the progress of one ribosome is obstructed 
by another. This picture is complicated by the fact that whilst there are only 20 common amino acids, 
in the model organism Saccharomyces Cerevisiae there are 41 species of tRNA, i.e., there is redundancy 
in the genetic code. For some amino acids there are multiple tRNAs; furthermore it is often that case 
that there is a highly abundant tRNA and a rare tRNA which code for the same amino acid. Sometimes 
a "slow codon", which will cause a long pause in elongation, is used when a fast codon corresponding 
to the same amino acid is available. The interest for biologists is therefore to understand how slow 
codons effect protein production rates and the use of ribosomes [B] . Translation reaches its completion in 
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Figure 1: Colour on-line, (a) Petri net components: places, transitions, arcs and tokens, (b) A simple 
Petri net. 

the termination stage where, through the binding of release factors, the ribosome disassociates from the 
mRNA and the amino acid chain is released ready for folding or further processing. The quantities of 
interest which we shall take from the model are the time interval between successive amino acid chains 
being completed, which we will call the protein production time, and the ribosome occupation density of 
the mRNA. 

In the next section we introduce Petri nets, and detail how these can be used to describe translation. 
Then in Sec. [3] we introduce the max plus algebra, giving a brief survey of the salient facts required in 
the rest of the paper. In Sec. |4]we detail the max-plus treatment of the Petri net describing translation. 
These results are then compared to numerical simulations: we first consider simple "designer mRNAs" 
where each codon corresponds to the same tRNA, before treating realistic sequences taken from the S. 
Cerevisiae genome. Finally we consider how introducing stochasticity into the model is likely to alter 
the results, and how the present work relates to previous models of translation (specifically the totally 
asymmetric simple exclusion process or TASEP). 

2 Timed Petri Nets 

Petri nets (attributed to C. A. Petri [7]) are a scheme where a sequence of discrete events is described on 
a network. We give a brief introduction here, but for a detailed description refer the reader to Ref. [5] 
and references therein. 

A Petri net is a directed graph which contains two types of nodes: places and transitions. These are 
shown diagrammatically in Fig. [IJa). Directed arcs connect places with transitions, but not places with 
places or transitions with transitions. Places can contain objects called tokens. A transition is said to 
be active when each of the upstream places it is connected to contains at least one token. Figure \V[b) 
shows an example Petri net. When a transition is active it can fire; upon firing one token is removed from 
all places upstream of the transition, and one token is added to all places downstream of the transition 
(clearly there is no implicit conservation of tokens since the number of upstream places need not equal 
the number of downstream places). Events unfold in a discrete manner. Petri nets are often used to 
model systems where events occur (transitions fire) given that a set of conditions are fulfilled (tokens are 
present) . 

A timed Petri net is an extension to this framework in which a waiting time is attached to each place. 
Whenever a token is put into a place, it will only contribute to the activation of transitions once a time 
t associated with that place has elapsed. This can be thought of as a timer on the token which is started 
as it enters a place. As soon as a transition becomes active it fires. In this way the Petri net describes 
discrete events occurring in continuous time as conditions arc fulfilled. Unlike some other models (for 
example Monte Carlo based simulations) events can occur simultaneously. 

We can describe the movement of ribosomes along the mRNA as a timed Petri net, and this is shown 
in Fig. [21a). A sequence of codons of length n is represented by two rows of places, with each pair (upper 
and lower) representing a codon. A token in the top row indicates that a ribosome is decoding that 
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Figure 2: Colour on-line, (a) A Petri net describing mRNA translation. Each pair (upper and lower) 
of places represent one codon. Tokens (red) in the bottom row represent vacant codons; in the top row 
they represent ribosomes reading that codon. The initial condition with no ribosomes on the mRNA is 
shown, (b) An example Petri net for a short mRNA of length n = 4 codon, with a token configuration 
representing ribosomes occupying codons 2 and 3. 



codon, and a token in the bottom row indicates a vacant codon. Having this two row structure ensures 
that there can only be one token in each place; i.e., one ribosome translating each codon. This of course 
requires that the initial state also satisfies this condition. A suitable initial state is having all top row 
places empty and all bottom row places containing one token; this corresponds to a mRNA initially free 
of ribosomes. 

The firing of a transition T; represents the movement of a ribosome from codon i to codon i + There 
are n + 1 transitions, labelled i = 0, . . . n, with the firing of the zeroth and nth transitions representing 
initiation and termination events respectively. At the leftmost end of the Petri net we have only an upper 
row place; this represents a continuous supply of ribosomes ready to begin translation. At the rightmost 
end of the Petri net there is only a lower row place, representing the cytoplasm which is always ready to 
accept ribosomes. These places are kept full with a token via an arc to and from the same transition. 

The time associated with each upper row place corresponds to the time it takes a ribosome to capture a 
tRNA of the species corresponding to that codon. Generally the arrival times of tRNAs can be represented 
by a set of independent Poisson processes, however a deterministic version of the model can be obtained by 
replacing these random times with their means. In the present paper we concentrate on the deterministic 
version; the consequences of introducing stochasticity are discussed in Sec. El We assume that the mean 
codon waiting time is proportional to the concentration of the species of tRNA corresponding to that 
codon. We denote the waiting time for the ith codon counting from the left Tj. We assume that when a 
ribosome moves forward the codon it was covering is immediately vacant, and the time it takes for the 
ribosome to physically move is negligible; thus we associate a time zero with each of the lower row of 
places. In reality a ribosome actually covers more than just the codon it is translating, but for simplicity 
we do not take this into account here. The first place in the chain - corresponding to initiation - has a 
time ri n associated with it. Biologically this corresponds to the average time interval between attempts 
for a ribosome to begin translation, and will depend on the availability of ribosomes and various initiation 
factors, as well as the presence of any secondary structures in the non-coding leader region of the mRNA. 
Since our focus here is on elongation, we assume that r m is constant in time, treating it as a control 
parameter. The final place in the chain has an associated time r out which is the time a ribosome waits at 
the end of the mRNA before releasing the completed protein and detaching from the mRNA. Biologically 
this could be linked to the availability of a number of release factors, and again we treat it as a control 
parameter. The model is therefore parametrised by the initiation and termination waiting times Tj n and 
r out , and the set of internal tRNA capture waiting times {t^I < i < n}. 

Figure [U(b) shows and example Petri net for a very (unrealistically) short mRNA of length n = 4 
codons, with ribosomes occupying codons 2 and 3. Transition T2 cannot fire because there is no token in 
its input place in the lower row (representing a vacancy). Transition T3 will fire as soon as the token in 
its upper row input place has been there for a time T3 . If the token in the place upstream of T2 (upper 
row place for codon 2) has been there for a time longer than ti when T3 fires, then T2 will also fire at 
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this time. Thus if there is a slow codon (large Tj), when the corresponding transition fires, it could result 
in several other transitions simultaneously firing. 



3 The max-plus semi-ring 

If each place in a Petri net has exactly one upstream and one downstream transition (as is the case in 
those depicted in Fig. [5]), then its dynamics can be described using max-plus algebra. Max-plus is an 
algebra over the semi-ring R max 6 MU-oo. In this section we give a brief overview of the theorems 
required in the rest of the paper, and refer the reader to Refs. [3l [8] for further details and proofs. 
In max-plus algebra the addition and multiplication operators © and © are defined as 

a (B b — max(a, b) and a ® b = a + 6, 

where a, b G tt mra , and we define e = — oc and e = which behave as zero and unity respectively. 
Max-plus algebra refers to the set 7?. max = (R max , ®, ©, e, e). The elements of 7vL max have much the same 
properties as in conventional algebra; for example 

a © b — b © a, 
a©(6ffic)=a(g>&ffia©c, 
a © e = e, 
a © e = a, etc. 

We define Max plus powers in the natural way 

a® x = a © a © ■ ■ ■ £g> a x times, 

and note that in conventional algebra this corresponds to x x a. We note that the <g> operator has an 
inverse which can be expressed as a negative power 



but there is no inverse of the © operator. We also use the following notation for sums and products over 
indices: 



cii = max{ai|l < i < n} and (^) > </,-. 



E 

i—l i—1 i—1 



Vectors and matrices are also defined, and we denote by [A]ij the i-jth component of the matrix A. 
The sum and product of matrices A s and B s 1^™ are then defined as 

[ A © B] ij = a ij © bij , 
/ 

[A © B]ij = £f) a ik © b k j , 
fe=i 

where a^- = [A]ij and bij — [B]ij. 

Matrices can have associated eigenvectors and eigenvalues; i.e. a vector u satisfying 

A © u = A © u, 

is an eigenvector of A associated with eigenvalue A. Note that eigenvectors are not unique, i.e. if u is an 
eigenvector then so is a © u with a an arbitrary finite number. In general a matrix can have multiple 
associated eigenvalues. 
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The properties of a max-plus matrix can often be determined by considering the associated directed 
weighted graph. The graph associated with matrix A, denoted G(A), consists of a set of nodes and 
directed weighted arcs, where if there is a matrix element [A]ij with value aij 7^ e, then there is an arc 
from node j i with weight (note that direction of the arc differs to that in conventional algebra and 
graph theory). A series of of one or more arcs between two nodes i and j is called a path from i — >• j, and 
if there is a path i — > i this is called a circuit and is denoted 7. A graph is said to be strongly connected 
if for any two different nodes there is a path between them, and a matrix A is said to be irreducible if 
the corresponding graph Q{A) is strongly connected. The circuit weight w 1 of a circuit 7 is defined as 
the sum of the weights of all arcs in that circuit, and the circuit length Z 7 as the number of arcs in the 
circuit. The mean circuit weight is defined w 7 — w 7 /l 7 . If the maximum mean circuit weight in a graph 
is A, then a circuit with a mean circuit weight equal to A is called a critical circuit. The critical graph 
corresponding to matrix A, denoted Q C {A), is defined as the sub-graph of G{A) containing only the nodes 
and arcs which are in the critical circuits. The cyclicity of a graph o~g is defined as the greatest common 
divisor of the lengths of all of the circuits in that graph, and the cyclicity of a matrix A, a a is equal to 
the cyclicity of the critical graph of A, G C {A). 

Theorem 1 Any irreducible matrix A G R^a™ possesses one and only one eigenvalue X, which is a finite 
number and is equal to the maximal mean circuit weight of circuits in the graph G{A). (See, for example, 
Theorem 2.9 in J3jj for proof.) 

Thus the eigenvalue of an irreducible matrix can be found by considering the corresponding graph. As 
noted above there is a whole continuum of eigenvectors associated with the eigenvalue of a max-plus 
matrix. Two eigenvectors x, y are co-linear if there exists a scalar a G K max such that y = a ® x, and 
such vectors can be projected onto the same object in a projective space [3]. A given eigenvalue of the 
matrix A can be associated with more that one linearly independent eigenvector, and the number of such 
eigenvectors can be found by considering the critical graph G C (A). 

Theorem 2 If the critical graph G C {A) of an irreducible matrix A has k maximal strongly connected sub- 
graphs (m.s.c.s.), then A has k linearly independent eigenvectors. (For proof see, for example. Theorem 

4 mm-) 

Linear equations such as 

x = (A®x)©b, (1) 
where A G M^ax an d x, b G Kmax: can often be solved using an object called the Kleene star, defined as 

00 

A*=Q)A® k . (2) 

k=0 

The existence of A* can be proven if the graph G(A) has only non-positive circuit weights [3]. 

Theorem 3 If G(A) has maximal mean circuit weight less than or equal to e — 0, then x = A* (Si h, is 
a unique solution to Eq. flj). (For proof see, for example, Theorem 2.10 in J3jj.) 

We shall see below that the problem studied in this paper involves sequences of vectors x(fc), k 6 N 
described by a recurrence equation 

x(fc) =A® x(fc - 1), 

or equivalently 

x(fc) =A® fe ®x(0), (3) 
for k > where A G R m a™ i s irreducible, and x(0) is some initialisation vector. 
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Theorem 4 The cyclicity theorem states that there is an integer Kq such that 

A ®(k+a A ) = A ®<x A ^ A ®k^ k > A - 0j 

where A ararf a a are the eigenvalue and cyclicity of the matrix A respectively. (For proof see or 
Theorem 3.9 in J$.) 

This implies some periodicity in the powers of A; the periodic behaviour is characterised by A and a a, 
and Kq is known as the transient time of A. By using this in Eq. we find 

x(k + a) = A^ k+ ^ ® x(0) 

= A 8<T ® A® k ® x(0) 

= \® a ®x{k). 

If matrix A has cyclicity a a = 1, then we have 

x(k + 1) = A (gi x(fc) 

= A®x(fc), fc>if , (4) 

i.e., x(fc) (or equivalently A® k ® x(0) for any x(0) £ t^) is an eigenvector of A for k > Kq. The effect 
of the initial condition has died out. 



4 Analysis of translation with the max-plus algebra 

A Petri net can be described using max plus algebra by defining the matrices Aq , A\ , . . . , Am which depend 
on the initial conditions such that [A m ]y is equal to the maximum of the holding times associated with 
the places between transitions j and i, which initially contain m tokens. The dynamics are then described 
by the vector x(/c), where the zth component of the vector is equal to the time at which the ith transition 
fires for the fcth time. This vector satisfies the recursion equation 

x(fc) = A (g> x(fc) © A x <g> x(fc - 1) 8 • • • © A M ® x(fc - M). 



For further details see [31 [T01 fTTj . 

In the case of the Petri net which describes mRNA translation (shown in Fig [U(a)), the structure and 
initial conditions are such that each place can contain either or 1 tokens, i.e., in this case M = 1 and 
the above equation reduces to 

x(k) = A ®x(fc) ®Ai (g>x(fc- 1). (5) 

In the Petri net we have n + 1 transitions labelled 0, . . . n, so x(fc) has components {xi\0 < i < n}. 
Likewise Aq and A\ are (n + 1) X (n + 1) matrices, with elements labelled by indices running from to 
1. These matrices are 



Aq 



I 



\ 



c 

Tl 

c 
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e 

T-2 



e \ 

c 
c 



( 



and A\ 



\ 



\ 





Tout / 



Written in component form Eq. ([5]) gives 

x (k) = max{r in + x (k - 1), xi(k - 1)}, 

Xi{k) = maxjri + Xi-i(k), Xi + \{k — 1)} for 1 < i < n, 

x n (k) = max{r„ + x n -i(k), r out + x n (k - 1)}. 



G 



If we think of the integer k as labelling the ribosomes (xi(k) gives the time at which the fcth ribosome 
leaves the ith codon) then these equations make sense conceptually since a ribosome will leave codon i 
either a time after it left codon i — 1, or at the time when the (k — l)th ribosome leaves codon i (i.e., 
time Xi+i(k — 1)), whichever is latest. 

Although the solution of the recursion equation §5§ is not straightforward, since there is no inverse to 
the © operator, by identifying the second term as a vector b = A\ <E> x(k — 1), we note that this equation 
is of the form of Eq. ((T|) , and so via Theorem [3] the solution is 

x(fc) = Aq <g> A\ <g> x(k - 1). 

From that theorem we know that Aq exists since the graph G(Aq) contains no circuits of positive weight; 
this can also be easily demonstrated by considering the first few powers of Aq [3J. For example for a 
system with four transitions (n = 3) we see 



e 
e 

Tl (g> T 2 



e 

t 2 ® r 3 



e e 
e e 
e e 



Aq 3 



c 
e 
e 

T"2 



T3 



i.e., the number of non-zero entries decreases as the power increases, and actually 

oo n 

A* - £D A® k — (X^ A® k 

fe=0 k=0 

The Petri net describing translation can therefore be described by the equation 

x(fe) = B<E>x(k- 1), 

where B = Aq ® A\ is given by 

/ T in e e 

Tin <g> Ti Ti £ 



(G) 



Tl ® T 2 



Tl T 2 



Tl 03 T 2 <X> T 3 Ti (3 T 2 



T3 



^2 



T 2 <K> T 3 







V 



Tl 



Tl 



T 2 



(7) 



'Tout / 



The eigenvalue and cyclicity of the matrix B can be found from the corresponding graph G(B). In 
Fig. [3Jwe show Q{B) for the example of n — 3. We note that the graph is strongly connected, so B 
is irreducible and has exactly one eigenvalue (Theorem [TJ . Due to non-zero elements on the diagonal 
of JB, there exist circuits of length one, so the cyclicity of the graph is og = 1. By inspection of the 
graph it is clear that one of the circuits of length one will always have a mean circuit weight ui 7 equal to 
the maximum mean circuit weight. That is to say, the maximum mean circuit weight, and therefore the 
eigenvalue of B will be equal to whichever of the waiting times is largest, i.e., 



A = max{ T in ,T out ,Ti 1 1 < i < n}. 
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It is also clear then that the critical graph of B must always contain a circuit of length one, so <jb = 1- 
Due to the cyclicity theorem (Theorem^]) , we conclude that for large k the vector x(fc) is an eigenvector 
of B. Regardless of the initial condition x(0), the Petri net describing translation always reaches a steady 
state described by 

x(fc) = A®x(fc- 1), (8) 

for k > Kq, where Kq is the transient time of A. The component of the eigenvector Xi(k) gives the time 
at which the ith transition fires for the fcth time. The nth component of this vector equation tells us 
that the transition corresponding to completion of an amino acid chain will fire at time intervals of A; 
thus the eigenvalue of A can be interpreted as the protein production time. As noted above, another way 
of interpreting the eigenvector is that Xi(k) is the time at which the kth token (ribosome) leaves the ith 
place (codon); the token enters the ith place at time Xi-\(k) and the maximum length of time which a 
ribosome can stay at a given codon is A. The proportion of time for which there is a token occupying the 
ith place in the steady state, i.e. the occupation density, is therefore given by 

pi = ~(xi -Xi-i) i = l,...n, (9) 

where {xi\0 < i < n} are components of any eigenvector of B, since Xi(k) — Xi-i(k) = Xi(l) — Xi~i(l) for 
all k, I > Kq. We note that pi are the components of a vector of length n. 

It now remains to evaluate the eigenvector. If we define u = x(fc — 1), then from Eq. ((SJ, x(fc) = A<g)u. 
We can then use this to evaluate the eigenvalue using Eq. ([S]) which when written as components gives 

A ® Uq = Ti n ®«o © «li (10) 

A ® Ui = A ® Tj ® Ui-i © Ui+i, i = l,...n — l, (11) 
A ® u„ = A ® r„ ® u n -i © r ou t ® u n . (12) 

Thus the eigenvectors takes a different form depending on the eigenvalue. We identify the following three 
cases: 

Case (i): A = r m 

We first consider the case where the parameters are such that the waiting time for initiation of translation 
is longer than that for tRNA capture and termination, i.e., r; n > r out ,Ti for 1 < i < n. We shall denote 
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this the entry limited regime, as it represents the situation where entry of the ribosomes into the system 
is the rate limiting process. Replacing A — > Tj n in Eqs. (|104I12|) gives 



T in ® Uq = Tin © U ffi Ml , 

Tin © U» = Tin © Tj © © U i+ l , 1 = ! 

Tin © «n = Tin © T n © T out © li n . 

From (|15|) . since T; n ^ T out , for consistency it must be the case that Ti, 
therefore 

U n = T n © U n -l. 



(13) 

.n-1, (14) 

(15) 

r„ © m„-i > T out © m„, and 
(16) 



From (fl"4"|) , taking 



n — 1 gives 



Tin ©"n-1 



T„_i © U„_ 2 © U n 



Since Ti„ ^ t„, in order for this to be consistent with Eq. (|16[) it must be the case that Ti n ©Tj ©w„_2 > 
meaning 



U n -2 



© «n-l 



Continuing to use Eq. (|14D in this way we can find expressions for each component of u in terms of the 
nth component. The eigenvector is given by 



®(-i) 



u 



T® ( '^lln \ 
7n © U n 



© M„ 



Using Eq. © we find = t^/t^. 
Case (ii): A = r out 

In the case where T ou t > Ti n ,Ti for 1 < i < n, the time it takes a ribosome to leave the mRNA is the 
limiting process, i.e., this is the exit limited regime. Following the same method as above we make the 
replacement A — » T out in Eqs. (fT2|) giving 

(17) 

...n-1, (18) 

(19) 



Tout © U = Tin © MO © "l, 

Tout © Ui — Tout © n © Ui-i ffi Ui+l , i 
Tout © U„ = Tout © T„ © U„_i © Tout © U n . 



This time we start with the first of the three equations, and note that since Ti n 7^ T ou t, it must be the 
case that u± > r- m © uo; therefore 

U\ = Tout © U . (20) 

We then consider Eq. (JTSJ); for i — 1 this gives 



Tout © Ul = Tout © Tl © U ffi 1t 2 . 
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In order for this to be consistent with Eq. (|20|) . it must be the case that U2 > T ou t ® t% ®uq, meaning 

U 2 = T ou t ® Ui 



We continue using (IT8|) to find the other components of the vector; the result is 

/ Uq \ 

Tout ® MO 

leading to occupation densities pi = 1. 
Case (iii): A = r 

The final case is the elongation limited regime, where one or more of the internal waiting times has a 
value t, with r > Ti n ,T ou t,{Tj =/= t}. We assume initially that there are two places with waiting times 
equal to t deep in the bulk of the mRNA, denoting the first the pth, and the second the gth. That is to 
say Tp, i~q — t where 1 < p < q < n. We again take Eqs. (|12[) and this time make the replacement A — > t, 
giving 

T ® Uq = Tin ® MO © «lj (21) 

t (g> m, = t ® Ti ® Uj-i © Mj+i, 1 < i < n, (22) 

T U„ = T (g) T„ ® © Tout <& U n - (23) 

We start with Eq. (|21L which since T; n 7^ r, gives 

Ul = T ® Uq. 

Moving onto Eq. (|2"2"|) and taking i = 1, if n 7^ r it must be the case that 

U2 = t ® Ui 

= T® 2 ® U0 . 

We can continue using Eq. (|22p in this way until we reach the pth waiting time, since t p = r. Taking 
i = p — 1 and then i = p in Eq. (I22p we find respectively 

r ® u p _i = u p 

t ®u p = t® 2 u p -i (B U p+ l. 

Since neither term on the right hand side of the second equation would contradict the first, all this can 
tell us is that 

up+i < t® 2 u p ^. (24) 
We now consider Eq. (f2"3"|); since T ou t 7^ t it must be the case that 



L n ' n 



<U n -i. 



(25) 



Using Eq. 



with i = n — 1 gives 



T ® U„-l = T ® T„_i ® M„_ 2 © U r , 
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If r n ^ r then in order not to contradict Eq. (|25[) . it must be the case that r eg> t„_i Cg> u„_2 > w„ 
Therefore 

u„-i = r„_i eg) it„_2- 



We can continue to iterate backwards using Eq. 
i = q and then i = q — 1 gives respectively 



until we reach the qth waiting time r q = r. Taking 



U q = T ® 

Again the second equation can only give an inequality 

Uq > T <g> T q -l ® M g _ 2 . (26) 

We are therefore left with two inequalities, Eqs. (|24l) and ([26| which cannot alone tell us about the pthe 
and the gth components of the eigenvector. We find that if we assume that one of these inequalities is 
actually an equality, then we can continue using Eq. (I22[) to find all of the components of the eigenvector 
consistently with all of the above equations. This leads to two eigenvectors depending on which inequality 
we set equal 
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where we recall that p < q and the pth and gth are the upstream and downstream most codons with 
waiting times equal to r. Any max-plus linear combination either of these two vectors is an eigenvector 
of B. If we consider the critical graph of B, if there are M non-adjacent codons with waiting times 
equal to r, then this graph will contain M m.s.c.s., and therefore (via Theorem [5]) there are M linearly 
independent eigenvectors. We note that the critical graphs in cases (i) and (ii) have only one m.s.c.s.. 
Each of the two eigenvectors above gives rise to a different occupation density profile using Eq. (J9j> , 

i.e., 

1 for 1 < i < q, 
Ti/r for q < i < n, 



and 



Pi = 



1 

Ti/r 



for 1 < i < p, 
for p < i < n, 



which corresponds to queues of ribosomes behind the qth and the pth codons respectively. In mRNAs 
with more than two isolated codons with waiting times equal to r, there will be solutions with queues 
behind each of these "slowest" codons. We will show in the next section that only one of the solutions is 
realised, and we reserve further discussion until Sec. 15.21 
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Figure 4: Colour on-line. Colour maps showing how (a) the mean ribosome density p and (b) protein 
production time P vary as a function of r; n and r out , for a uniform mRNA where all codons have the 
same tRNA capture time r, = f . All times are in units of f . Plot (c) shows the results for P again, but 
also includes in blue the boundaries between the three regimes as determined from the max-plus algebra. 



5 Simulation Results 

In this section we compare the results from the max-plus algebra detailed above with direct numerical 
simulation of the Petri net shown in Fig. [U]a). We conduct the simulations by considering an initial state 
where all bottom row places contain a token and all top row places are empty (the first and last places in 
the chain both contain one token). We then move forward in time from one event (firing of a transition) 
to the next, updating token positions at each step. We allow the system to reach a steady state, and 
then examine the protein production time P, and the codon occupation density of each place pi, as well 
as the mean density, defined as 

1 " 

P=-E«- (27) 

i=l 

We first consider a very simple mRNA where each codon is assumed to be identical, i.e., it codes for 
the same tRNA, and the waiting time for each is the same. The internal hopping rates are chosen to be 
Tj = f for 1 < i < n, and f is used as the unit of time. We then measure P, pi and p for different values 
of T in and r out . 

In Sec. I5.2l we consider real mRNA sequences from the S. Cerevisiae genome. The waiting time for a 
particular tRNA species is assumed to be inversely proportional to the concentration of those molecules 
found in a typical cell. We estimate this from the gene copy number of each tRNA [12]. We choose the 
internal hopping rates such that the average of the waiting times is equal to f , and again use this as the 
unit of time. 

5.1 Homogeneous mRNAs 

Figure [T] shows colour maps for the protein production time and mean density as a function of Ti n and 
T ou t for a uniform mRNA of length n = 500, as generated from numerical simulations. There are three 
regimes, with the dynamics depending on the relative magnitudes of f , T m , and r ou t- These correspond 
to the three cases found from the max-plus algebra, i.e., the initiation, termination and tRNA capture 
limited regimes. We summarise the results in Table [TJ Figure [5] shows how P and p vary as a function of 
(a) Ti n for fixed T ou t, and (b) T ou t for fixed T{ a . The numerical and analytic results match exactly. 
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(i) Initiation limited 


(ii) Termination limited 


(iii) tRNA capture limited 




(Tin > Tout , f ) 


(Tout > r in ,f) 


(f > Tin, Tout) 


Protein production time P 


Tin 


Tout 


f 


Occupation density pi 


T/Tin 


1 


1 


Mean density p 


T/Tin 


1 


1 



Table 1: Details of each regime as given by the max-plus algebra, for a uniform mRNA, where the tRNA 
capture waiting times are Tj = f for 1 < i < n. These are found to exactly match the numerical simulation 
results. 



(a) (b) 




0.5 1 1.5 2 0.5 1 1.5 2 

Tin T ou t 

Figure 5: Colour on-line. Comparison of numerical and analytic results for a uniform mRNA. Plots 
showing (a) the mean density p and (b) the protein production time P at different values of T m and r ou t . 
Squares show numerical results for r out = 1.5f, and crosses those for r out — 0.5f at different values of 
Ti n . Red lines show the analytic results from the max-plus algebra. All times are shown in units of f . 

5.2 Heterogeneous mRNAs 

We now consider mRNA sequences from genes YJL136C and YDR382W in S. Cerevisiae, which we here 
on denote mRNA A and mRNA B, respectively. For a given mRNA of length n codons, each of the 
waiting times < i < n} takes a value chosen from the set of 41 times {si\l < i < 41} corresponding 
to each of the 41 tRNA/codon species, and (1/41) J2t=i Si ~ r - We take f as the unit of time; biological 
experiments estimate that f k 0.1 s -1 pQ. Figure [5] shows the waiting times at each codon for each of 
the sequences. 

Figures [7] and \E\ show simulation results for P and p for each mRNA respectively. We again observe 
three regimes depending on the relative magnitude of Ti n , T out , and the largest of the Ti which we denote 
t. Again the simulation results exactly match those from the max-plus algebra, which are summarised 
in Tabled] 

In two of the three regimes the occupation density differs from codon to codon. In Fig. [9] we show the 
density profile for each mRNA when parameters are chosen such that the system is in the entry limited 
phase. We note that for most codons i, pi < 0.5; in other models (i.e., the TASEP which we discuss in 
Sec. [7]) the equivalent regime is usually called the low density phase. In Fig. [10] we show the density 
profiles for the tRNA capture limited phase. For codons located upstream of the most upstream codon 
with the longest waiting time (t p — t) the density pi ~ 1 for i < p; i.e., the leftmost instance of the 
slowest codon causes queueing of ribosomes. In the case of mRNA B the first "slowest codon" is at i = 1, 
so no queue is observed. Although from the max-plus algebra there are other possible solutions, only 
one is realised in the simulations. Instances of the slowest codon appearing further downstream do not 
give rise to queues, since the rate at which ribosomes arrive at those codons is too low. In terms of the 
max plus algebra, the different possible eigenvectors are reachable by different initial conditions [15] . The 
initial condition where are the codons are vacant corresponds to the initial vector x(0) where xq(0) = 
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Figure 6: Colour on-line. Waiting times for each codon on mRNAs A (top) and B (bottom). Red dots 
indicate codons which have the largest waiting time, which for mRNA A is r = 0.7696f and for mRNA 
B is r = 0.6155f. Note that each waiting time Ti has a value from the set {si\l < i < 41}. Whilst 
(V 41 ) Eiii s * = f > in general (1/n) £" =1 n f. 





(i) Initiation limited 

(Tin > Tout , t) 


(ii) Termination limited 

(Tout > Tin, r) 


(iii) tRNA capture limited 

(t > Tin, Tout) 


Protein production time P 


Tin 


Tout 


T 


Occupation density pi 


Ti/Vin 


1 


1 for 1 < i < p, 
ti/t for p < i < n 


Mean density p 


1 " 

— 

«Tin f— ' 
i— 1 


1 


n 

- + > Ti 

n nr z — ' 



Table 2: Details of each phase as given by the max-plus algebra, for a non-uniform mRNA. Here r is the 
longest internal waiting time (r = max{rj|l < i < n}), and p is the position of the leftmost codon with 
waiting time equal to r. 
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Figure 7: Colour on-line. Comparison of numerical and analytic results for mRNA A showing values of 
p and P for different values of Ti n and r out . In (a) the three regimes are clearly visible. In (b) and (c) 
squares show numerical results for r out = 1.125f and crosses those for r out = 0.375f. Lines show the 
analytic results from the max-plus algebra. For this mRNA r = 0.7696f. 
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Figure 8: Colour on-line. Comparison of numerical and analytic results for mRNA B showing values 
of p and P for different values of T ln and r out . In (a) the three regimes are clearly visible. In (b) and 
(c) squares show numerical results for r out = 0.75f and crosses those for r out = 0.25f. Lines show the 
analytic results from the max-plus algebra. For this mRNA r = 0.6155f. 
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Figure 9: Colour on-line. Occupation density profiles for each mRNA when parameters are chosen such 
that the system is in the entry limited phase. Parameters for mRNA A are T; n = 1.3125f, r out = 
0.0375f , r = 0.7969f , and for mRNA B r in = 0.875f , r out = 0.025f, r = 0.6155f. Bars show numerical 
results and crosses show results from the max-plus algebra. Red dots indicate the locations of the slowest 
codons. 



and Xj(0) = e for i ^ 0; i.e. we start counting the elapsed time from at the first initiation event. Only 
the eigenvectors corresponding to a queue behind the most upstream slowest codon are reachable from 
this initial condition. 

The other solution for the density profile of mRNA B can be accessed by making a change to the 
waiting times during the simulation. If we first allow the system to reach steady state, and then increase 
the waiting time on the slow codon at position i — 104 by some small amount 8t, then a queue will begin 
to form behind this codon. If we reverse the change after the system reaches steady state, then the queue 
behind codon 104 will persist (data not shown). 

6 Introduction of stochasticity 

In the above work we have examined a deterministic version of the model, where rather than choosing 
waiting times from a Poisson distribution, we instead use the mean of the distribution. Simulation of a 
Petri net with waiting times chosen from a distribution is straightforward. The description using max-plus 
algebra is less so. At its essence the problem consists of a product of irreducible max-plus matrices such 
that x(fc + 1) = A(k) ®x(fc) for k > 0, where {A(k)\k £ N} is an independent and identically distributed 
(i.d.d.) sequence of random matrices. Some work on such sequences can be found in the literature [3lll4j. 
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Figure 10: Colour on-line. Occupation density profiles for each mRNA when parameters are chosen such 
that the system is in the tRNA capture limited phase. Parameters for mRNA A are Tj n = 0.0375f , r out = 
0.0375f , t = 0.7969f , and for mRNA B r in = 0.025f, r out = 0.025f, r = 0.6155f. Bars show numerical 
results and crosses show results from the max-plus algebra. Red dots indicate the locations of the slowest 
codons. 



and it can be shown that for an i.d.d. sequence of matrices which have fixed support and are irreducible, 
a max-plus Lyapunov exponent exists, and this is equal to the asymptotic growth rate. However, there 
is currently no method for calculating this quantity, which in the present model would be related to the 
mean protein production time. For this reason, here we only briefly discuss how simulation results for a 
stochastic Petri net compare with the deterministic case. 

We consider here a uniform mRNA of length n = 500 codons, where waiting times for initiation, 
termination and tRNA capture are chosen from exponential distributions with means T m , r out and f 
respectively. We perform simulations where we allow the system to reach steady state before measuring 
the time averaged protein production time and codon occupation density, which we again denote P and p 
respectively. As shown in Figs. Illf a - ) and (d) we again see three regimes depending on which waiting time 
is the largest (r in and T out compared to f). We note by comparing with Fig. 0]that in the stochastic case 
the boundary for the tRNA capture limited regime is at r m , r out ~ If , compared to r- m , r out = f in the 
deterministic case. Also in that regime, the density is significantly reduced compared to the deterministic 
case (p w 0.75 compared to p — 1) and the protein production time is increased (P ~ 2.7? compared to 
P = f). This is as expected since allowing the internal waiting times to vary will lead to gaps between 
ribosomes. The fact that the difference between the stochastic and deterministic models in this regime 
is so large can be attributed to the fact that the rate limiting process involves choosing waiting times for 
each of the n = 500 places representing the codons. In the other regimes the effect of the stochasticity 
is less severe, as the limiting process involves only one place. Deep within the initiation or termination 
limited regimes p and P are almost the same as in the deterministic model. The general effect of stochastic 
waiting times in these regimes is that both density and protein production time increase. If we consider 
ribosomes occupying two consecutive codons, and if the waiting time drawn for the leading ribosome is 
long, it could hold up the ribosome behind it; if the time drawn for the leading ribosome is short it is 
unlikely to effect the ribosome behind it. Thus we expect ribosomes will in general move more slowly 
along the mRNA; as we shall see below, often the maximum of two exponentially distributed random 
numbers is what determines the overall behaviour. 

For non-uniform mRNAs the effect of stochasticity in the initiation and termination limited regimes 
are similar: both p and P are increased compared to the deterministic model. In the tRNA capture regime 
the situation is more complicated, and features such as multiple queues are observed. Discussion of such 
phenomena is beyond the scope of the present paper. Here we only briefly consider one consequence 
of stochasticity which has played an important role in other translation models [15l [16]. It has been 
observed in a widely used alternative model for translation (the TASEP, as discussed in section [7]) that 
several slow codons in close proximity have a more dramatic effect on the protein production rate from 
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Figure 11: Colour on-line. Simulation results for a stochastic Petri net representing a uniform mRNA 
of length n — 500 codons, where each codon corresponds to the same tRNA species. Plots (a) and (d) 
show colour plots of the mean occupation density and protein production time interval respectively for 
different values of r in and r out . The other plots show line graphs at cross sections through the colour 
maps as indicated by the dashed blue lines. Dashed lines in (b), (c), (e) and (f) show results for the 
corresponding deterministic Petri net. All times are given in units of the mean tRNA capture time f. 
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Figure 12: Colour on-line. Schematic diagram showing a cluster of two slow codons labelled A and B. 
Ribosomes translate from left to right. All other codons have much shorter waiting times compared to 
A and B. 



an mRNA than slow codons in isolation. In the deterministic Petri net the token release time interval 
from a place corresponding to a slow codon (a slow place) which is in isolation is the same as that from 
the second of a pair of adjacent identical slow places. Once steady state has been reached the "timer" 
on the token in the first slow place will end at the same time as on the second, so a token will always 
enter the second place as soon as the previous token leaves - the time interval between tokens leaving the 
second place is exactly the waiting time for that place. The deterministic Petri net does not predict any 
effects due to clusters of slow codons. For a simple toy mRNA sequence it is possible to quantitatively 
predict how a cluster affects the protein production rate in the stochastic version of the Petri net model. 

We consider an mRNA with very fast initiation and termination containing only very fast codons, 
except for a single slow codon somewhere in the bulk. We focus on the release time interval of tokens 
from the place corresponding to the slow codon. If the slow place has mean waiting time r s i ow , and we 
assume that all of the other waiting times are so fast that this place is refilled with a new token practically 
immediately after the previous one leaves, then the mean release time interval for tokens from this place 
is equal to r s i ow . Now consider the same system, but with a pair of slow codons labelled A and B - see 
Fig. HU In the Petri net we denote the places corresponding to the slow codons place A and B. We denote 
the waiting times drawn from the (identical) distributions for each place Ty, t^, . . . , and rf 1 , t% , 
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We set t = when the system is in the steady state and a token enters place B (since the other codons 
are very fast, a token also enters place A at this time). The first token leaves place B at time t = Tj 8 ; the 
next token leaves B at t = + maxjr-j 4 , Tj 8 }, and the next at t — + max{r^, tv 8 } + maxjrj 4 , r-f}. 
That is, the ith token leaves place B at 

i-l 

* = if +^max{r/,rf }. 

.7=1 

The ith interval A< between two consecutive tokens leaving B is therefore 

A i = rf -r^+max^.r^}. 
Taking the average, and noting that the distribution for each slow place is the same, gives 

(A. 1 ) = (max{r^ 1 ,r i B 1 }), 

i.e., the mean release time interval is the mean of the maximum of two numbers drawn consecutively 
from the same exponential distribution. This can easily be showrQ to be 3r s i ow /2, where r s i ow is the 
mean waiting time of the slow codons. This agrees with our earlier assertion that, in general, stochastic 
waiting times leads to slower movement of ribosomes along the mRNA. 

7 Relationship to other translation models: TASEP 

An alternative framework for describing mRNA translation, is the totally asymmetric simple exclusion 
process, or TASEP jTTJ [TS1 [IS]. In that model ribosomes are represented by particles which hop in one 
direction along a ID lattice of sites, which represent the codons of the mRNA. The current model shows 
many similarities to the TASEP, but the dynamics unfold in an essentially different way. 

The dynamics of the TASEP are most often simulated via Monte Carlo methods using a random- 
sequential update rule, or using continuous time Monte Carlo |19j . In the former case, time is discretised 
and in each time step lattice sites are chosen at random with uniform probability. If the chosen site 
contains a particle and the next site is vacant the particle is advanced with some probability p] particles 
are injected from the leftmost site and removed from the rightmost site with probabilities a and /3, 
respectively. For a lattice with n sites this is repeated n times in each time step. In the continuous time 
version, exactly one particle is moved in each time step but the particle is chosen based on its probability 
of hopping, and the duration of the time step drawn from the corresponding Poisson distribution. Both 
methods are equivalent in that they are the realisation of the usual master equation in continuous time; 
this last point means that multiplying all probabilities by a common factor leads only to a rescaling of 
time. Under some conditions the random sequential TASEP can be solved exactly |18j . and also a mean 
field treatment which ignore spatial correlations is often used. 

The major difference between this and our timed Petri net picture is that in the TASEP only one 
particle can move at a certain instant in time. In the Petri net, as soon as any token move becomes 
possible, it is executed; i.e., events can happen simultaneously. Since in translation ribosomes actually 
take a finite time to move from one codon to the next, i.e., they are able to move at the same time, 
we propose that the Petri net more closely models the microscopic dynamics of this system. Results 
from the Petri net model can most easily be compared with the TASEP by identifying aSt — 1/Tj n , 
/36t = l/r out and, in the case of a uniform mRNA, p5t = 1/f, where St is the time step. Although the 

x If two times t\ and T2 are random numbers where the probability that T{ has a value x is given by a Poisson distribution 
with mean t s i ow , i.e., , P(t; = x) = e~ a! /" r siow ! then the probability that m = max{Ti,T2} has a value x is given by 
P(m = x) = P(ri = x) ■ P(t2 < x) + P(t2 = x) ■ P(ri < x). The mean value of m is then just the expectation value of 
P(m = x). 
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Table 3: Comparison of the results from the Petri net and TASEP models for a uniform mRNA. The 
TASEP results are from a mean field treatment which has been shown to very closely match random- 
sequential Monte Carlo simulations for large systems [18] . 

uniform TASEP produces a very similar phase diagram to that of Fig 01 the subtle differences in the 
dynamics do lead to some important differences in the macroscopic behaviour. The analytic results from 
the max-plus treatment for the Petri net model and from a mean field approximation for the TASEP are 
summarised in [3] As well as the density and protein production time differing in the two models, the 
boundaries of the regimes are also different. Specifically the boundaries for the tRNA capture limited 
regime (usually called the maximal current phase in the TASEP literature) are at T; n ,T ou t = 2f instead 
of f in the deterministic Petri net. We note that the phase boundaries in the TASEP are the same as 
those in the stochastic Petri net model discussed in the previous section; we could therefore ascribe this 
to the stochasticity. However in the tRNA capture limited regime we have p = 1/2 and P = 4f in the 
TASEP, compared to p = 1 and P = f in the deterministic Petri net and p ~ 3/4 and P « 2.7 f in the 
stochastic Petri net. The dynamics of the TASEP are still different to those of the stochastic Petri net. 

A TASEP with an ordered-sequential update has been studied by some authors [20] [HJ HI] ; and is 
particularly applicable to, for example, traffic flow models, movement of molecular motors, and - as we 
study here - mRNA translation. This update rule corresponds more closely to the Petri net model. Sites 
are taken in turn from right to left and updated according to the hopping probabilities. However, this has 
remained a less favourable update rule since analytic approaches run into several difficulties, particularly 
in the physical interpretation of hopping probabilities. Unlike the random-sequential case, for the ordered 
update time is not continuous, i.e., scaling the probabilities does not lead to a simple scaling of time (for 
further discussion see [H]). The advantage of the Petri net picture is therefore clear: we can include the 
fact that multiple events can happen simultaneously, whilst retaining (at least in the deterministic case) 
an analytically soluble mathematical framework, i.e., algebra on the max-plus semi-ring. 

Another difference between the current model and the TASEP is in the biological interpretation of the 
waiting time of the ribosome on each codon. In the random-sequential TASEP the hopping probability 
from a site is chosen based on the abundance of tRNAs which corresponds to that codon. A recent 
study |23] shows that this is in fact the wrong interpretation of the hopping probability. Ciandrini et 
al. have developed an extension to the TASEP where ribosomes take two internal states. They identify 
two times: the waiting time for a ribosome to capture the correct tRNA T cap t U re, and the time it takes 
for the ribosome to physically move from one site to the next r movo . The former of these depends on 
the availability of tRNAs whilst the latter does not, and they argue that the capture of the tRNA can 
occur independently of whether or not there is a vacancy to the right of the ribosome. A commonly held 
misconception in applying the TASEP to translation is that the model represents the limit r movo — > 0; 
Ciandrini's model shows that it is the opposite T cap turc — > limit which recovers the original TASEP 
model; importantly, this limit does not describe the biologically relevant regime (r movc <C T cap t U ro)- In 
the timed Petri net the place waiting times correspond to T cap t U re, hence we explicitly operate in the 
biologically relevant T movo — > limit. 
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8 Concluding Remarks 



We have presented here a new model of ribosome dynamics during mRNA translation. Algebra on the 
max-plus semi ring lends its self to describing systems in which discrete events occur depending on the 
fulfilment of conditions, and has previously be used to study, for example, distributed software systems, 
automated manufacturing or industrial control systems. We have shown here that biological systems 
represent another area where such methods can be applied. The Petri net is a useful tool for visually 
representing discrete event systems, and the wealth of previous work on describing Petri nets using max- 
plus has allowed us to quickly develop a framework for predicting the dynamical behaviour of elongating 
ribosomes given only the codon sequence of an mRNA. 

The analytic treatment we have presented is applicable to a deterministic version of the model, and 
we have discussed the modifications to the system when ribosome waiting times are chosen from a 
distribution. The present work therefore represents the first step in a new direction for protein synthesis 
modelling. Whilst numerical simulation of a stochastic Petri net is straightforward, there are currently 
few methods or algorithms available for the study of sequences of i.d.d. max-plus matrices. Therefore 
this work also represents a new source of motivation in this area. 

There is also scope for introducing features which bring the model closer to the biology. For example 
including a finite time for the physical movement of the ribosomes |23| , taking into account the fact that 
ribosomes cover more than one codon at a time [23, 25 , or allowing for a finite pool of ribosomes or other 
resources [26l [27l [28 ] . 

In summary we have presented a new model of translation which uses max-plus algebra to solve a 
Petri net description of the system. The microscopic dynamics are more realistic than in other translation 
models, and calculation of quantities such as density profiles and protein production rates is straightfor- 
ward. Unlike other translation models the max-plus algebra also allows for exact analytic solutions for 
inhomogeneous mRNA sequences. 
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